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As part of our ongoing lattice study of the electric polarizabilities of hadrons using the background 
field approach, we use reweighting to examine the effect of the field on the sea quarks. As with 
other reweighting studies, the chief difficulty lies in the construction of a stochastic estimate of the 
ratio of the fermion determinants. In contrast to the case of reweighting in the quark mass, these 
estimators converge extremely slowly, and are resistant to common variance-reduction techniques 
such as low-mode subtraction. However, it is possible to construct an alternate estimator, taking 
advantage of the fact that we are interested in only perturbatively small fields; this estimator is 
susceptible to a variance-reduction technique based on a hopping parameter expansion. 
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1. Introduction 

Recently there has been a great deal of interest in first-principles lattice computations of hadron 
polarizabilities. The most interesting such quantity is the electric polarizability of the neutron, since 
it is difficult to access experimentally due to the lack of free neutrons; the best measurements have 
been obtained by neutron- lead [1] and neutron-deuteron [2] scattering. A basic computation of 
this quantity is not too hard to do using the background-field method [3]; the difficulties come in 
the approach to the physical point that can be compared with experiment. Effective field theories 
predict very strong dependence on m q near the chiral limit, making the neutron polarizability a 
good probe of whether or not chiral behavior is accurately captured by lattice simulations; doing 
so requires lighter quark masses and potentially expensive chiral actions [4]. Similarly, we must 
address finite- volume effects [5] and finally extrapolate to the continuum. 

The most difficult effect to accommodate, however, is the effect of the electric field on the 
sea quarks. In principle, this could be done from the ground up by including the background field 
(see Sec. 3) in gauge generation itself. However, these ensembles will necessarily be uncorrelated. 
To determine the polarizability we examine the mass shift in the neutron (or other hadron) when 
a small background electric field is applied; since the zero-field and finite-field propagators are 
strongly correlated, the error on this mass shift can be much smaller than the error on the masses 
themselves. Comparing two uncorrelated ensembles destroys this correlation. What we would like 
is to generate two (or more) correlated ensembles with different values of E. This can be achieved 
via reweighting. 

The most difficult aspect of a reweighting calculation is stochastic estimation of the weight 
factors; the present work is chiefly concerned with this problem. For reweighting in the quark 
masses, a straightforward stochastic estimator coupled with several standard variance reduction 
techniques is generally successful in obtaining good estimates with a reasonable amount of com- 
puter power. These techniques, however, are not useful for reweighting in the background field. We 
will discuss possible reasons for those failures, which will likely be illuminating for other groups 
performing reweighting computations. We then present an alternative approach in which we con- 
struct a stochastic estimator for the derivatives of the weight factor with respect to the background 
field which is susceptible to an improvement technique based on a hopping parameter expansion. 

The only prior work known to us on sea contributions to polarizability was done with a purely 
perturbative method [6]. Instead of applying a uniform background field of specified strength 
throughout the lattice and then computing hadron propagators, the author first expands the path in- 
tegral in powers of E and then computes the needed diagrams on the lattice using current insertions. 
In this method the disconnected and connected insertions (corresponding to sea and valence con- 
tributions) appear quite naturally. The author applies this method in a mixed-action formulation, 
computing domain-wall valence propagators on Asqtad dynamical configurations generated by the 
MILC collaboration [7] and using stochastic estimators to compute the disconnected diagrams. 

2. Lattice simulation details 

We apply the following methods to a series of gauge ensembles with two flavors of nHYP- 
smeared clover fermions [8] with m % ~ 330 MeV and a standard Symanzik-improved gauge action 
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with /3 = 7.1, giving a = 0.1255 fm (determined by the Sommer scale ro) [9]. The ensembles 
have volumes 24 3 x 48, 32 x 24 2 x 48, and 48 x 24 2 x 48; these elongated lattices were originally 
generated for a scattering study [10], but we reuse them here, since the elongation in the x\ direction 
is a convenient probe for finite- volume effects associated with Dirichlet boundary conditions. Each 
ensemble has 300 minimally-autocorrelated configurations. 

3. The background field method 

The most commonly-used approach for computing hadron polarizabilities is the background 
field method, allowing the extraction of the polarizability from spectroscopic measurements. The 
effect of a uniform background field on the mass of a hadron can be parametrized as: 



where jU is the magnetic moment, p the electric dipole moment, a is the electric polarizability, j8 
is the magnetic polarizability, and the ellipsis includes various higher-order terms as well as spin 
polarizabilities [11, 12]. While we are mostly interested in the neutron electric polarizability, the 
reweighting approach we will use to probe the effects of the sea is not specific to any particular 
type of hadron. The basic approach, then, to extract the polarizability of some hadron is to measure 
its mass both in the presence and the absence of a perturbatively small electric field (chosen small 
enough that higher-order effects are small) and examine the mass shift. 

To apply a background electromagnetic field to the lattice, one can simply apply a U(l) phase 
on the gauge links on top of the dynamical SU(3) gauge configurations, making the transformation 



We choose here to apply a constant electric field in the x\ direction. This can be done in any 
suitable gauge; we choose A4 = iEx\. The factor of i arises from the Wick rotation to Euclidean 
time; thus, on the lattice, a real electric field gives a real factor e^ Xi , while an imaginary field gives 
pure phase factor. Provided that the magnitude of the field is small enough, it should be possible 
to perform a lattice calculation for an imaginary field and analytically continue the results to the 
real axis; preliminary studies have confirmed that both methods result in consistent results for the 
polarizability [12]. It is convenient to use an imaginary field because the links remain unitary and 
the Dirac operator remains 75-Hermitian. Thus, we apply the electric field by the transformation 
U4 — > e~ ,r]Xl where 77 = a 2 qE. Note that this value depends on the charge of the quark flavor in 
question, so to compute the neutron correlator we will need gauge links for two values of 77 . 

While in the limit of infinite statistics there is no order-77 shift in the hadron mass, one may 
appear from statistical fluctuations in a real calculation. To eliminate this source of error, we com- 
pute hadron correlators for fields +f]x\ and —f]x\ (which should, in the infinite-statistics ensemble 
average, be identical) and take their geometric mean to get the two-point function in the presence 
of the electric field G(t, 77) = yj T\)G— (?, 77). To extract the polarizability, we then fit it along 
with the zero-field correlator G(t,0) to the form 




(3.1) 



(3.2) 



G(t,r])=Ae- {MN+8l]1) 



(3.3) 
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to get the mass shift 8. The two correlators must be fit jointly since they are drawn from the 
same gauge configurations and their fluctuations are (strongly) correlated; this correlation greatly 
reduces the statistical error in 8. 

3.1 Boundary conditions 

Periodic boundary conditions are preferable for many lattice observables. However, they po- 
tentially create substantial problems for measurements using background fields. For the gauge we 
have chosen and for real electric fields, there will be an unavoidable and large discontinuity in the 
electric field at the lattice boundary; for imaginary fields, this discontinuity can be avoided only by 
choosing particular values of E [13]. However, such fields are out of the perturbative regime that 
we want to examine. For other choices of the gauge, other problems manifest themselves, such as 
an electric scalar potential that is not single- valued, leading to manifestly nonphysical scenarios in- 
volving quark lines winding around the torus in the direction of the electric field. It is not clear how 
one should address these effects. We can avoid them by applying Dirichlet boundary conditions in 
the x\ and X4 directions. Dirichlet boundary conditions create their own problems, of course; we 
can now no longer use certain improved nucleon sources, such as Coulomb-wall, and can no longer 
project completely onto a zero-momentum state. However, these effects can be considered to be 
finite-size effects which go away in the infinite volume limit, and are manageable in the analysis; 
we thus use Dirichlet boundary conditions [12]. 

3.2 Choosing a field strength 

We must choose a value of the parameter v\ to use for the background field. (Since there 
are two quark flavors, two T]'s are required.) Choosing a value which is too large means that we 
leave the perturbative regime and begin to probe 6(E A ) effects; choosing a value which is too 
small means that we may encounter issues with numerical precision, either with the accuracy of 
inverters or (in the extreme case) machine precision. Fig. 1 shows the response of the neutral pion 
correlator to the background field at different temporal separations as a function of the d quark 
77 for an inverter precision of 10~ 13 ; the onset of nonperturbative behavior is clear. We note that 
when a similar - study was done with an inverter precision of 10~ 9 , the expected quadratic scaling 
behavior broke down at the smallest T7 's. This illustrates the need to choose an 77 small enough to 
avoid higher-order effects when probing polarizabilities, and the need to use a sufficient inverter 
accuracy to avoid numerical artifacts. Within the large flat region in Fig. 1, we are free to choose 
whatever value of rj makes the reweighting process perform best. 



4. Reweighting 

Reweighting is a technique for extracting physics based on a different action than the one 
used in Monte Carlo ensemble generation; essentially, it allows for post hoc modification of the 
parameters in the action. In the standard quantum Monte Carlo, we would like to do a path integral 
of the form 

J[dU]0e- s ° (4 1} 



J[dU]e- s ° 
where So is the QCD action. 
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Figure 1: Dependence of (log G ^ ^ 
different correlator times. This quantity is roughly equivalent to the shift in the effective mass divided by 
rj 2 , and should be constant in the range where rj creates a purely quadratic effect. The breakdown of the 
quadratic behavior is evident at large 77 . 

By generating a Monte Carlo ensemble with each configuration weighted by e~ So , the path 
integral reduces to the familiar 

Eg _ L# 
£1 ~ N 

Suppose, however, that we wish to use this ensemble (weighted by e~ So ) to learn about the behavior 
in the presence of the background field, given by the action Sjj . Then the Monte Carlo average gives 



(4.2) 



(4.3) 



£ e -(S,,-So) 

where e~( 5 i~ s o) is a "weight factor" which may be interpreted as the relative prominence of a 
particular configuration in the target and source distribution. If there is little overlap between 
the source and target ensembles, reweighting will fail as the weight factors fluctuate wildly (over 
many orders of magnitude); in general, reweighting always comes with a decrease in statistical 
power, as described in [14]. This may not always be immediately evident and it is possible to wind 
up underestimating statistical errors [15]. Reweighting is used to reweight in the quark mass to 
approach the chiral limit without incurring the expense of simulating at those light quark masses 
directly [15], to perform simulations at small but nonzero chemical potential, to couple sea quarks 
to dynamical photon fields to probe isospin breaking [16], and (most similarly to this project) been 
used to investigate the intrinsic strangeness of the nucleon [17]. In this last case, the authors sought 



e ^ by measuring Mn at slightly different values of m s and examining the difference. 



This is difficult if the errors in are uncorrected, but by reweighting in the strange quark mass 
to generate correlated ensembles they were able to measure errors on the difference in more 
accurately than Mjy itself. We intend to do a very similar thing, except with v\ (which can be 
thought of as either reweighting in the background field or the quark electric charge) instead of in 
the strange quark mass. 



4.1 Estimating the weight factor 

The weight factor w,- is given by e^ Sr >^ s °^ = det^ = det _1 M M n l , where M and M r] are 
the fermion matrices corresponding to the actions Sq and S^. This determinant is impractical to 
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compute exactly and must be estimated stochastically. The standard stochastic estimator for the 
inverse determinant of some matrix Q. is [15] 



While the inverse determinant of MqM^ is real, using this operator with the above estimator 
will produce complex results. We thus instead apply the stochastic estimator to an operator with the 



same determinant but which is Hermitian: Q. = yM^~ MqMqM^ . The square root can be taken 
by a rational function approximation as done in rational hybrid Monte Carlo, etc. This is somewhat 
computationally expensive, as it requires inverting a matrix which itself requires the computation of 
inversions; however, as that matrix is in principle perturbatively close to 1 its inversion should not 
require too many iterations. As an alternative, the authors of [1 8] have argued that in such cases one 
should simply discard the imaginary part of the estimator, saving all the work associated with the 
rational function approximation. In the analysis that follows, we estimate the inverse determinant 



While this stochastic estimator may be (and generally is) quite noisy, this noise does not in- 
troduce bias because the average over the noise vectors t, commutes with the gauge average [15]; 
in principle, if the ensemble is large enough, one need only use a single stochastic estimate per 
configuration, although with a gauge ensemble of finite size it is often profitable to work harder 
than this to improve the stochastic estimator to reduce its fluctuations. The simplest way to do 
this is to average multiple stochastic estimates, but there are other techniques that can result in a 
greater reduction in the stochastic noise for a given budget of computational power. We note that 
the "signal" we are trying to extract from these estimates is the true fluctuation of the weight factor 
from one configuration to the next, and thus we may adopt the rough criterion for reduction of the 
stochastic noise that CJg aU ge > Cnoise- This condition is not necessary, but it should be sufficient, for 
a successful application of reweighting. (Note, however, that if the true values of the weight factor 
are indeed very similar, it is possible for this test to give an overly-pessimistic description of the 
quality of the stochastic estimates.) 

5. Improving the stochastic estimate 

There is a fundamental problem, however: this estimator is tremendously noisy in our case. 
As one might expect, the distribution of stochastic estimates e~^= is log-normal; when the 

width is large enough that the skew of the distribution is apparent, sampling the long tail becomes 
very difficult. For small values of tj , the average estimate of the weight factor is very close to unity, 
while for larger values the long-tail sampling problem becomes tremendously difficult; see Fig. 2. 1 
We have confirmed that the estimator does produce the correct result for the determinant on a 4 4 
lattice where that determinant can be computed exactly [19], but it required 10 5 -10 6 noise vectors, 
something clearly unaffordable for production-size lattices! 

'While it is of course possible that the true value of the weight factor on the configuration shown here is close to 
unity, tests on other configurations give the same result. 




(4.4) 





M Mn where w,- = det Q. 
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Figure 2: Histograms of the stochastic estimator of the weight factor for three different values of 

rj, along with the mean and its standard error, on a single 24 3 x 48 configuration. Even with a large number 
of noises, the average is indistinguishable from unity; for large values of tj, the average is dominated by the 
few points in the right tail of the distribution, which is roughly log-normal. 

5.1 Low-mode subtraction 

A commonly-used technique in reweighting in the quark mass is low-mode subtraction. By 
projecting out the low modes of the operator in question and computing their inverse determinant 
exactly, we stand to improve the stochastic estimator in two ways. First, we accelerate the inver- 
sions required to make the estimates by reducing the conditioning number of the operator being 
inverted, allowing for the computation of a larger number of estimates with the same amount of 
computer time. Second, if an appreciable fraction of the fluctuations of the estimator come from 
the low modes, we may eliminate those fluctuations by treating the low sector exactly (by simply 
multiplying together the eigenvalues). This is offset, of course, by the additional cost of computing 
the eigensystem. Can this technique be profitably applied to reweighting in the sea quark charges? 

Computation of the eigensystem is more difficult in our case. In the mass-reweighting case, 
the eigenvectors for all values of fc are the same, so the eigenvectors of Q. are the same as the eigen- 
vectors of the Dirac operator. On the other hand, we must do a separate calculation to compute the 
eigensystem of D. = J M^~ 1 MqMoM^ 1 , which is quite computationally expensive. Furthermore, 
the eigenvectors differ for different values of TJ, causing problems with the determinant breakup 
method (see Sec. 5.2). 

Since our matrix £2 is in principle close to the identity, we computed eigensystems for both 
high and low sectors to test the method, generated an ensemble of a few hundred noise vectors, and 
used those to estimate det -1 Q. varying the size of the extremal sectors treated exactly. Specifically, 
if P is a projector onto the space spanned by N extremal modes of Q. with eigenvalues A,-, and 
P = 1 — P, we may write 

(detar^^n^^^)- (5 - 1} 

To our surprise, there is essentially no effect whatsoever on the size of the stochastic fluctuations, 
as shown in Fig. 3. This is strikingly different from the behavior for mass reweighting [15]. 

In mass reweighting, most of the signal comes from the low modes, as shown empirically 
in [17]; in fact, almost the entirety of their signal comes from the low modes, and while the high 
modes are included for correctness their contribution seems to be only (a small amount of) noise. 
The authors of [20] show this analytically, and furthermore argue that in the case of mass reweight- 
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Q.= 



1 MqMqM T i 1 for r\ = 0.00018, showing pairing between high and low modes. 



ing most of the noise comes from the low modes as well; the estimator is protected from large 
fluctuations from the high modes. Thus, by removing the low modes from the stochastic estimator, 
it is possible to also remove most of the noise. This is only guaranteed in the case of reweighting 
in the quark mass, however; in our case, it fails. Since changing the quark mass does not change 
the eigenvectors, the low modes of the Dirac operator itself are also eigenvectors of the matrix Q. 
for quark mass reweighting; this is not true for us, and the extremal modes of Q. may not be related 
to the low modes of the Dirac operator. 

For reweighting in the sea quark charge, the extremal modes of Q. do not contribute much to 
the determinant. This is because they are very nearly paired; for each high mode with eigenvalue 
A, there is a low mode with an eigenvalue close to nearly cancelling their contribution to the 
determinant; a representative case is shown in Fig. 3. Studies on 4 4 and 6 4 lattices confirm that 
most of the signal (the difference of the determinant from unity) in our case comes from the bulk 
"interior" modes of Q. 

5.2 Determinant breakup 

Another commonly-used technique for improving the stochastic estimator is to divide the 
reweighting up into many steps and compute an independent estimate of each one; this technique 
is often referred to as "determinant breakup" [15] or "determinant factorization" [14, 21]. For 
reweighting in the quark charge, for instance, one might imagine estimating the weight factor to 
reweight from T7 = to T7 = TJi, then the weight factor from T7 = t?i to tj = T?2, and so forth until the 
desired value is reached. However, for reweighting in tj , this technique also fails to decrease the 
stochastic noise in the estimator compared to simply using more noise vectors. The authors of [14] 
find that for reweighting in a substantial shift in m;, it is more efficient to use more subintervals in 
the determinant breakup (which they call "steps") than repetitions of the entire procedure ("hits") 
for a given total number of inversions. However, this improvement is limited, as shown in [18], 
and, as predicted in [14], there is no more to be gained after A^teps increases past a certain point. We 
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suspect that this behavior occurs because determinant breakup obtains its benefits by converting a 
large reweighting step into a sequence of small ones; once the steps are already sufficiently small, 
no further benefit is gained by splitting them up further compared to simply increasing A^its- There 
is no benefit to applying this technique to reweighting in 77, however, and this explains why; we are 
already reweighting by a perturbatively small interval in order to make a valid measurement of the 
polarizability. 



6. The perturbative estimator 

Since it is not feasible to perform a direct stochastic estimation of the determinant, and since 
neither of the most common improvement techniques that are successful for mass reweighting work 
for reweighting in the quark charges, we turn to an alternative estimator. We are interested only in 
the weight factors for perturbatively small values of 77, so we can expand about 77 = 0: 



w(r\) = 1 + 7] 



dw 
drj 



+ 



1 , d 2 w 



17=0 



2" 



(6.1) 



17=0 



The expansion must be taken to second order in 77, since we are interested in effects of order E 2 . 
The task then becomes to construct estimators for and . Both terms are needed, 

dr \ 77=0 W TJ=0 

since the linear term in the weight factor can couple with a linear dependence on 77 to give a 
quadratic effect, or the quadratic term can give a second-order effect on its own. Given estimates 
of these derivatives, we can evaluate the above at any sufficiently-small 77 to produce a reweighted 
ensemble on which to apply the valence calculation. 

, we rewrite detM^ as a Grassmann 



To construct such an estimator for |^ 



integral: 



77=0 



d detM,, 
drj detMo 



77=0 



detM 



d_ 
drj 



I 



dydy 



dM 



dM. 



detMTr ( ^AT 1 
\dr] 



(6.2) 



This is an expression for 



d detM 



, but we want 



-detM 



drj , uui worn detMo drj 

77 = 0, the determinant in front cancels the and we get 



77=0 



When the above is evaluated at 



d detM, 



d 77 detMo 



= Tr 



?7=0 




M 



-1 



(6.3) 



?7=0 



The trace involved here must still be evaluated stochastically. 



9 



Sea Polarizability via Reweighting 



Walter Freeman 



0.03 



e- 0.02 



o.oi 




0.00 




10" 



10" 



10" 
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The second derivative term proceeds similarly. Computing -^-j detM, we get 



dr] 1 



dr\ : 



- detM 



dr] 1 
d 



d\\fd\\fe 



-if/My 



d\\fd\\f 



dM 
- ¥ — ¥ , c 



dydy -y 



d 2 M 



+ j dv/dvj j Y^W 



-if/My 



Evaluating these integrals gives: 



- detM = detM 
drj z 

As above, this gives 

d 2 detM 



Tr 



d 2 M. 
drf 



r M~ 



dM 
Tr— AT 
drj 



(dM 



drj 2 detM,, 



Tr (BMq 1 ) + Tr (AM l AM 1 ) - (Tr (AM 1 ) ) ' 



where A 



dM 
dr] 



andfi: 



rj=0 



r)=0 

d 2 M 
drj 1 



77=0 



(6.4) 
(6.5) 



(6.6) 



(6.7) 



Now we must construct stochastic estimators for these three traces, as well as the trace in the 
first derivative. We use the standard stochastic estimator Tr ^ = (E^ note that two stochastic 
estimates of the first derivative term can be used to construct an estimate of (Tr(AM -1 )) 2 . Tests 
on a 4 4 lattice confirm that this estimator for the first derivative is correct; see Fig. 4. However, it 
is just as noisy as the direct estimator of the determinant; the data in Fig. 4 required 5 x 10 6 noise 
vectors to produce! Clearly this new estimator on its own has not gained anything. 

6.1 Hopping parameter expansion improvement 

However, this estimator for the trace is susceptible to an improvement technique. If other op- 
erators G\ can be identified such that the stochastic fluctuations in t, f and t, ' are correlated, 
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then we can reduce the overall fluctuations by writing 



Tr0=(tf <A+£Tr^. 



(6.8) 



Obviously, for this to be useful, the ^-'s must themselves have exactly-computable traces. 

For the operators needed here, a set of improvement operators can be gotten by performing 
a hopping parameter expansion of M~ l , with each term in the expansion acting as one &[; this 
guarantees that the fluctuations in the estimator will be correlated. Specifically, the fluctuations are 
correlated with the magnitude of the off-diagonal elements of ^M~ l ; the terms in the hopping 
parameter expansion approximate the largest of these close to the diagonal [22]. For the Wilson- 
clover fermions considered here, we may write M as 1 — k(@ + C), where the hopping term 3> is 
given by 



where L^ v (n) is a sum of the imaginary part of the plaquettes in the /iv plane that include site n. 
For the first derivative term, the technique of 6.8 gives 



where G[ = 1 k'(S ! + C)'. The stochastic component of this can be computed directly as 

before. The difficult part is to compute Tr k" $M {3> + C) n quickly. In the end we wish to apply this 
technique to nHYP-smeared Wilson-clover fermions, but it is useful to consider the unimproved 
Wilson action first, for which C = 0. 

The trace of ^ & l can be evaluated most simply by writing £F as the sum of eight hopping 
terms. ^ can be calculated analytically, and consists itself of two hopping terms in the ±£4 
direction; the rest are zero. We break and @" into separate hopping terms and expand the 
product. All terms which give a nonzero contribution to the trace are products of hops which form 
closed paths. However, not all such paths contribute; some closed paths (such as those that double 
back on themselves) have zero Dirac trace and do not contribute. The most efficient way to compute 
the trace is thus to expand as a sum of products of hopping terms which have definite spatial 

and Dirac structure, each corresponding to the product of gauge links along a path of a particular 
shape, and pick out only those paths which are closed. We then compute the Dirac trace of each 
term. The computationally-intensive step, summing the appropriate products of gauge links over 
the spatial volume, must only be done for terms with nonzero Dirac trace. 



@mn = £ [(1 - 7n)U^(m)8 m+(l . n + (1 + 7/z)^(«)5 mj „+/i] 



(6.9) 



and the clover term C is given by 



C, m = CsWg^Oyiv^jUV^) 



(6.10) 
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The clover term can be dealt with by considering a generalization of the procedure used for 
the pure Wilson action. Another way to describe that procedure is as a separation of the Q operator 
into eight pieces (hops in ±x,y,z,t), each with a definite Dirac and spatial structure. The need 
to separate out terms with a definite spatial structure is readily apparent since in taking the trace 
we are only concerned with closed paths. In doing so we got a definite Dirac structure that does 
not depend on the gauge links and can be factored out of the sum over sites along the way. The 
same principle applies to the clover term even though the spatial structure is trivial: separate C into 
six pieces each with definite Dirac structure. These are the individual terms o^yL^y in Eq. 6.10, 
allowing us to factorize the trace into Dirac and SU(3) pieces. 

We must also compute as it appears in This can be done analytically for "plain" 
clover fermions, but for nHYP-smeared clover fermions, the smearing process complicates com- 
putation of |f and While this can in principle still be done analytically, it is simpler and not 
expensive to do it numerically. Our procedure for evaluating Tr K n ^ (J2> + C) n is thus as follows: 

1. Write 3l and |f (the latter calculated numerically) as the sum of eight separate hopping 
terms. Similarly, break C and |^ into six terms, each with the Dirac structure of a^ v . 

2. Expand ^§±2(#+C)», giving 14" +1 terms, each with a definite Dirac and spatial structure. 

3. Compute spatial part: The majority of these terms correspond to paths that are not closed, 
and can be discarded right away. 

4. Each term now can be factorized into a SU(3) part (which depends on the gauge links, and 
must be summed over sites), and a Dirac part (which is the same for each lattice site and can 
be factored out). 

5. Compute Dirac part: the majority of terms will have zero Dirac trace and can be discarded. 

6. Compute SU(3) part: For the remaining terms (roughly one in 500, regardless of order, after 
the first few), do the hard work of computing products of links and L^ v 's over all sites. 



The second-derivative traces can be evaluated in a very similar way. For the Tr [BM 1 ) term, we 
must additionally compute B = — y 2 ' , but this is not difficult to do numerically. 

This expansion can be in principle carried out to arbitrarily high order. The evaluation of the 
stochastic estimates of additional <^/'s is not difficult, as the cost is dominated by the inversion 
which only must be performed once for any number of orders. The limiting factor is the evaluation 
of TrK'^S}' , as the cost increases exponentially with i. However, it is possible to examine the 
reduction in stochastic noise without computing the exact traces. We find that the degree to which 
this improvement procedure reduces the stochastic noise, especially after the first few orders, is 
strongly dependent on the value of m n ; this is not surprising in light of the use of a hopping 
parameter expansion. This behavior is shown for several different values of m % in Fig. 5. Because 
the improvement slows down after the first few orders and the cost increases substantially, we 
compute the exact traces only up to seventh order. 
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7. Applicability to the polarizability and future plans 

While the improvement technique carried out to seventh order yields a substantial reduction 
in the variance of the stochastic estimator, this estimator is still quite expensive to run; we aim 
for around 10 3 stochastic sources per configuration. These estimates are almost complete on the 
24 3 x 48 ensemble. Once they are complete, we will use the estimates of |^ and |^ to evaluate 

vv(tj) = jg 1 ^ at particular values of rj (corresponding to the charges of the up and down sea 
quarks), chosen small enough that they lie within the perturbative regime of Fig. 1 but large enough 
to avoid any subsequent numerical precision issues, and then use the resulting weight factors to 
generate a reweighted ensemble in the presence of the background field. The remainder of the 
polarizability computation proceeds as in Ref. [3], but using the reweighted ensemble to compute 
the nonzero-field neutron correlator. It is likely that the influence of the charged sea is greater at 
lighter quark masses. We have such an ensemble available (with m K pa 200 MeV) and intend to 
repeat the calculation on it, although the variance reduction from the hopping-parameter expansion 
will not be as strong due to the lower m K . Finally, once the calculation detailed above is complete, 
if the variance in the stochastic estimator of the weight factors leads to a large increase in the 
overall error, we are considering using deflation to accelerate the thousands of inversions required 
to further improve the estimator. 
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